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ABSTRACT 

We establish constraints on the mass and abundance of black holes in the Galactic halo 
by determining their impact on globular clusters which are conventionally considered 
to be little evolved. Using detailed Monte Carlo simulations, and simple evolutionary 
models, we argue that black holes with masses Mhh ~ (1 — 3) X 10 6 M Q can comprise 
no more than a fraction fhh ~ 0.17 of the total halo density at Galactocentric radius 
R &t 8 kpc. This bound arises from requiring stability of the cluster mass function. A 
more restrictive bound may be derived if we demand that the probability of destruction 
of any given, low mass (M c « (2.5 — 7.5) x 10 4 M Q ) globular cluster not exceed 50%; 
this bound is fbh ^S 0.025 — 0.5 at R rs 8 kpc. This constraint improves those based on 
disk heating and dynamical friction arguments as well as current lensing results. At 
smaller radius, the constraint on fbh strengthens, while, at larger radius, an increased 
fraction of black holes is allowed. 
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1 INTRODUCTION 

What is the form and structure of dark matter in galactic 
halos? A variety of both baryonic and non-baryonic can- 
didates exist (see Carr 1994 for a review of baryonic dark 
matter candidates) but there are relatively few constraints 
, so the question remains. 

One longstanding suggestion is that of Lacey & Ostriker 
(1985) who proposed that halo dark matter consists of mas- 
sive black holes with Mbh ~2x 1O 6 M0. In so doing, they 
cast a solution to two problems: 1) what is the composition 
of the dark matter; 2) and what is the mechanism which 
heats the Galactic disk? Their calculation showed that a 
steady flux of 2 x 1O 6 M0 black holes passing through the 
disk would heat the disk in the manner required to explain 

1 /2 

the velocity dispersion-age relation <r* oc tj for disk stars 
(Wielen 1977). 

Although subsequent observational and theoretical 
work suggests that an explanation of disk heating does not 
require massive black holes (Carlberg et al 1985; Stromgren 
1987; Gomez et al 1990; Lacey 1991)- indeed, analysis of the 
disk heating problem is ongoing (e.g. Sellwood, Nelson & 
Tremaine 1998)- one can, in any case, view the disk heating 
argument as a disk heating constraint. The constraint can 
be developed by generalizing the Lacey & Ostriker model 
to Mth > 2 x 1O 6 M0 with a less-than-unity fraction of halo 
mass in black holes fbh < 1 (e.g. Carr, Bond & Arnett 1984; 



Wasserman & Salpeter 1994). Then, since the energy input 
to the disk AE oc M^ h for a single black hole, any combina- 
tion fbh Mbh ~2x 1O 6 M0 produces the same net heating of 
the disk. Therefore fbhM b h >2x 10 6 M Q overheats the disk 
and is definitely not allowed. The generalization fbh < 1 is 
desirable, given the variety of dark matter candidates, the 
results of microlensing surveys (e.g. Alcock et al. 1997) and 
the fact that not all dark matter need be baryonic given the 
bound from primordial nucelosynthesis (Pagel 1997). 

Are there other constraints on the mass and abundance 
of black holes in the Galactic halo? In some sense, halo black 
holes are suprisingly difficult to detect, given that there is 
considerable observational evidence for black holes of similar 
mass (10 6 M© < Mbh J$ 10 9 M Q ) in the centers of galaxies 
(e.g. Kormendy & Richstone 1995). Conversely, they have 
been surprisingly difficult to constrain or rule out. Lacey & 
Ostriker themselves remarked that the accretion luminosity 
of such objects may be too high to have escaped detection; 
however, no definitive constraint has been recorded (Carr, 
Bond & Arnett 1984; Carr 1994). Hut & Rees (1992) argued 
that dynamical friction would drag ~ 100 of these objects 
into the Galactic center in a Hubble time, leading to coales- 
cence and production of a central object much larger than 
allowed by observational constraints {Mbh ~ 2x 1O 6 M0; 
Genzel et al 1997). However, Xu & Ostriker (1994) tested 
this argument with detailed N-body simulations and found 
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that typically only one black hole would remain in the Galac- 
tic center due to three-body encounters. 

Constraints from gravitational lensing are compara- 
tively weak at this time. For the values of Mbh consid- 
ered in this paper, the large size of the Einstein ring gives 
event durations orders of magnitude too long for the present 
Galactic microlensing surveys. Lensing of quasars (Canizares 
1982; Kassiola, Kovner, & Blandford 1991) restricts Qbh = 
{pbh) / pc to be less than about 10%, which is greater than 
the estimated mass in dark haloes-an upper bound for the 
scenario we are considering. Several observing plans have 
been proposed to detect massive black holes. Turner and 
Umemura (1997) argue that the Sloan Digital Sky Survery 
and the Hubble Deep Field can place constraints on Qbh by 
looking for extremely large amplifications of O and G stars 
at cosmological distances. Also, Turner, Wardle, & Schnei- 
der (1990) propose that Mbh ~ 1O 6 M black holes are de- 
tectable through arcsecond size lensing of objects in M31 
and the Galactic Center. 

Wielen (1985,1988) first pointed out that globular clus- 
ters also constrain the properties of massive black holes in 
the Galactic halo because of their susceptibility to external 
heating and tidal disruption. Later, Moore (1993) applied 
the same arguments to a set of low mass globular cluster in 
the halo. These arguments have been re-examined in more 
detail by Klessen & Burkert (1995) and Arras & Wasserman 
(1998), who first included cluster evolution due to black hole 
heating in examining the constraints. Our goal in this and 
subsequent work will be to re-examine constraints on fbh 
and Mbh imposed by globular clusters using Fokker-Planck 
calculations of cluster evolution which include the effects of 
encounters with massive black holes. 

Both Wielen (1985,1988) and Arras & Wasserman 
(1998) delineate two mass regimes for black holes: the low 
Mbh regime, in which individual collisions perturb a clus- 
ter only weakly, but where many such collisions produce 
a steady, diffusive energy input; and the high Mbh regime, 
where a single encounter can destroy the cluster. In the low- 
mass limit, one can obtain limits on the product fbhMbh', in 
the high-mass limit, one can obtain limits only on fbh for 
Mbh > Mhigh- Applying these arguments to Moore's (1993) 
cluster sample, Arras & Wasserman (1998) concluded that 
fbhMbh S; 1O 3 M in the low-mass regime and fbh ;$ 0.3 in 
the range 1O 6 M < M bh < 1O 7 M . 

Although Arras & Wasserman (1998) included evolu- 
tion due to black hole heating, they did not consider the 
influence of internal relaxation and post-collapse evolution 
in their calculations. Thus they pointed out the need for im- 
proved calculations to derive the most robust constraints on 
fbh and Mbh using the globular cluster argument. 

To address this issue in the present work, we combine 
the statistical framework developed by Arras & Wasserman 
(1998) with the Fokker-Planck evolutionary calculations em- 
ployed by Murali & Weinberg (1997a-c). Our calculations in- 
clude two-body relaxation, post core-collapse evolution and 
tidal shocking by massive black holes. Using a Monte Carlo 
approach, we focus on the high-mass regime and directly 
compute the probabilities for strong collisions between glob- 
ular clusters and massive black holes. 

In order to translate these calculations of the evolu- 
tion of individual clusters into bounds on fbh, we need to 
consider the implications of our results for the evolution 



of a population of clusters. To accomplish this, we adopt 
two different points of view. The first, and more restric- 
tive, viewpoint is that global studies of the evolution of 
the globular cluster population are consistent with observa- 
tions assuming no black holes at all (e.g. Murali & Weinberg 
1997b), so that including black holes should have only a min- 
imal effect. Following this approach, we find that to ensure 
50% survival probability for globular clusters with masses 
M c w (2.5 - 7.5) x 10 4 M at R = 8 kpc, the fraction of 
the halo in black holes with masses Mbh <; (1 — 3) x 10 6 M 
must be fbh ,$ 0.025 — 0.05. This limit on fbh is between one 
and two orders of magnitude stronger than the disk heating 
constraint at this Mbh, and would imply that black holes 
with Mbh ~ 10 6 M are not a candidate for baryonic dark 
matter in galactic halos. 

Our second approach, which turns out to be systemat- 
ically less restrictive, examines the evolution and stability 
of the globular cluster population in the context of a sim- 
ple model. As we shall see, the two principal effects of per- 
turbations by black holes are to destroy clusters outright, 
and to cause surviving clusters to lose mass. We model this 
simply via a partial differential equation that includes mass 
loss via an advection term, as well as destruction. In order 
to obtain results valid for black hole masses ~ 10 6 M , we 
must extend the tidal approximation employed throughout 
most of this paper, so that non-destructive encounters at 
impact parameters inside the tidal radius of a cluster may 
be taken into account. (The extension of our calculations 
to include such penetrating encounters, as well as more de- 
tailed evolutionary models, will be treated by us elsewhere.) 
In the context of these models, we find that requiring the 
observed mass distribution of Galactic globular clusters to 
be stable over ~ 10 Gyr requires fbh J$ 0.17 at R — 8 kpc 
for M bh > 2 x 1O 6 M . This bound, although not as tight 
as our more restrictive (and more qualitative) limit on fbh, 
still implies that massive black holes cannot be the primary 
constituent of the halo dark matter. 

The plan of the paper is as follows. We summarize the 
framework for determining constraints in §2. In §3 we deter- 
mine which globular clusters may be considered relatively 
unevolved over the age of the Galaxy in the absence of bom- 
bardment by black holes, and then find the probability that 
they are destroyed for given values of Mbh and fbh- The re- 
sults are extended to different cluster and black hole param- 
eters using scaling arguments. We then discuss the proper- 
ties of clusters which are not destroyed outright, but instead 
undergo many non-destructive collisions over the age of the 
galaxy. Lastly, in §4 we discuss how our results constrain 
Mbh and f bh . 



2 FRAMEWORK 

To re-examine constraints on fbh and Mbh set by globu- 
lar clusters, we incorporate collisions and encounters with 
black holes into multi-mass Fokker-Planck calculations of 
cluster evolution, which include two-body relaxation and 
phenomenological binary heating of the core. Our code de- 
scends from that of Chernoff & Weinberg (1990). In prac- 
tice, we take clusters on circular orbits at R — 16 kpc: this 
minimizes the effect of relaxation and allows us to neglect 
Galactic tidal heating while subjecting clusters to the ap- 
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proximate black hole-flux crossing the disk. In addition, this 
radius corresponds to the spatial region containing many of 
the clusters in the sample used by Moore (1993). See Table 
1 below for a list of input parameters for the calculation. 

Given that clusters evolve and some will vanish through 
evolution (e.g. most recently Murali & Weinberg 1997a-c; 
Gnedin & Ostriker 1997; Vesperini 1997), it is important at 
the outset to establish which clusters to use in setting con- 
straints. The lifetime for given cluster orbit scales roughly 
with internal dynamical time tdyn and cluster mass M c as 
tuf e oc td yn M c for quasistatic evolution. Since td yn scales 
with Galactocentric position due to tidal limitation, low- 
mass clusters in the inner Galaxy are the first to vanish. 
In general, for any particular orbit, there is a minimum 
mass cluster which survives to the present-day. Clusters with 
evaporation timescales less than a Hubble time tn do not 
provide straightforward constraints on fbh and Mw,: clus- 
ters currently at or below the minimum mass might have 
had considerably larger initial mass. 

It is important to note that the predictions from evo- 
lutionary calculations appear quantitatively consistent with 
observations. The main uncertainty in the Fokker-Planck 
calculations is the core heating term: calculations predict 
that clusters have high central densities in the post core- 
collapse phase, while it is unclear precisely how this relates 
to observations (e.g. Drukier, Fahlman & Richer 1992). Nev- 
ertheless, the predicted death rates are not wildly incon- 
sistent with observations and differences of opinion arise 
mainly over the importance of evolution in clusters at the 
peak of the luminosity function, M c ~ 1O 5 M0 (Murali & 
Weinberg 1997b; Gnedin 1997; Harris et al 1998; Kundu et 
al 1998). 

With this in mind, we adopt a two-step approach to 
investigating constraints on black hole masses: 1) we first 
determine cluster initial conditions which do not strongly 
evolve in smooth halos in a Hubble time; 2) we then immerse 
these clusters in halos with black holes to investigate the 
possible constraints. While cluster survival is most likely for 
large M c , significant perturbation by black holes is most 
likely for small M c . We consider cluster masses near the 
lowest M c that can survive for 10 Gyr when fbh = 0. 

2.1 Dynamics of encounters 

We specify isotropic distributions of perturbing black holes 
using the fbh-Mbh parameterization: therefore the local 
number density of black holes nbh(R) = fbhphaio(R) / 'Mbh- 
This implies the relative velocity distribution for encounters 
given in equation (16) of Arras & Wasserman (1998). 

Encounters between clusters and black holes are pre- 
dominantly impulsive given the Galactic rotation velocity 
V c ~ 220 kms -1 . For a cluster at R ~ 16kpc, the charac- 
teristic internal velocity Vint ~ 5kms _1 . For a cluster on a 
circular orbit at V c and a random perturber drawn from a 
non-rotating, isothermal halo, the typical relative encounter 
velocity V re i ~ V c . For the black hole masses and abun- 
dances considered below, the influence of impacts at ~ 10 
tidal radii, rt, is small. Even at this distance, the timescales 
rt/vint >> lOrt/Vc so the probability of a non-impulsive 
encounter is very small. 

To perform the simulations described below, we incor- 
porate individual collisions between black holes and globular 



clusters into Fokker-Planck calculations using the impulse 
approximation. To determine the effect of each encounter, 
we calculate the second-order change in the distribution 
function (e.g. Murali & Weinberg 1997a). The method is 
similar to procedures used by Murali & Weinberg (1997a) 
and Gnedin & Ostriker (1997) in studies of cluster evolution 
in the Galactic tidal field. However, in the appendix, we show 
that there is a small error in the previous treatments. Our 
current treatment remedies this. 

This approach represents a linearization of the full col- 
lision problem. In complete generality, the collision problem 
requires simultaneous solution of the coupled, collisionless 
Boltzmann-Poisson equations. Linearization imposes a lim- 
ited range of validity. In the appendix, we show that this 
approach is valid for dM/M < 0.15 in a single encounter. 
We terminate any calculation where dM/M > 0.15. 

We adopt the Fokker-Planck approach, rather than, say, 
N-body simulations because we can study evolution due 
to both internal and external effects and because we need 
a computationally feasible method to conduct the Monte 
Carlo simulations described below. While N-body simula- 
tions permit fully non-linear calculations of strong collisions, 
it is difficult to include two-body relaxation, core heating 
which leads to post core-collapse evolution and the effect 
of weak encounters in which only small mass loss occurs. 
N-body simulations are also much too expensive to use in 
Monte Carlo simulations. 



2.2 High- mass limit 

Our analysis focuses on the high-mass limit for halo black 
holes (e.g. Bahcall, Hut & Tremaine 1980; Wielen 1985; Ar- 
ras & Wasserman 1998). The limiting mass is defined as the 
mass for which a single, tidal encounter at the typical rel- 
ative velocity can destroy a cluster. For completeness, we 
sketch the definition of the high-mass limit following the 
detailed discussion given by Arras & Wasserman (1998). 

Let us assume that cluster destruction occurs when a 
strong collision unbinds a fraction dM of the total mass. 
As shown in Appendix A, the fractional mass loss in the 
impulsive tidal limit 



\dM/M\ =f = KM2 h /V?*b* 



(1) 



where b and V re i are the impact parameter and relative ve- 
locity of the collision. The constant 



K 



KG 2 {r 2 } 



(2) 



The quantity (r 2 ) is the mean square radius of the cluster, 
a is its one-dimensional central velocity dispersion and k is 
a dimensionless constant depending only on its structure; 
see Arras & Wasserman 1998, equations (36) and (37). Note 
that, from the virial theorem, a 2 oc r// 1 and that, by tidal 
limitation in an isothermal sphere, rt oc R ' , so that K oc 
R 2 , where R is the Galactocentric radius. 

Rewriting this, we may define the 'destructive radius': 



bd= ( r mL 



1/4 



(3) 



where fd is the fractional mass loss leading to destruction. 
Given a black hole of mass Mbh moving at velocity V re i 
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with respect to the cluster, an encounter with any impact 
parameter b < bd leads to fractional mass loss f > fd from 
the cluster. Of course, since the black hole can travel with a 
range of relative velocities in relation to the cluster, we may 
invert this relation to define the range in V re i which leads 
to destructive encounters within rt: 



V re l< 



KMu 



1/2 



V, 



rel.max 



(4) 



In general, cluster properties depend on time, so that bd and 

V Te l,max will tOO. 

We define the limiting mass Mbh = Mhigh to be the 
black hole mass for which fd = 0.15 in an encounter at 
b = n with V re i = 2.5{V re i), where (V re i) w 1.4714. The 
factor of 2.5 is introduced to ensure that the probability of 
non-destructive collisions inside rt is very small. (The factor 
of 2.5 is probably larger than it needs to be: we estimate 
a probability of about 10~ 5 for a nondestructive collision 
inside rt when Mw, = Mhi g h with this choice.) Note that 
(bd/r t ) 2 = M b h/Mhi g h for this value of Vrei, and, more gen- 
erally, (bd/nf = 2.5{V re i)M bh /V re iM Mgh ). 



2.3 Monte Carlo enumeration of collision 
probabilities 

To proceed, we introduce a framework for calculating the 
probability that a cluster experiences an encounter with 
fractional mass loss fd in a Hubble time. For the black 
hole background specified above, the rate of encounters with 
/ > fd is 

= 2k I dV rel dbbV re ,n bh (R)F(V rel ,t) = r d , (5) 

/o Jo 

where F(V re i,t) is the relative velocity distribution of the 
black hole population with respect to the instantaneous mo- 
tion of the cluster. For a non-evolving cluster on a circular 
orbit, Arras & Wasserman (1998) show that 

1/2 

(6) 



dNd 
dt 



V, = -:/>;,;,!/.')( j 



Implicit in equation (6) are a number of noteworthy scal- 
ings: Nd is actually independent of Mbh > Mhigh and M c , 
and Nd oc _R _1 given fbh- In general, the probability that a 
cluster suffers an encounter with / > fd in a Hubble time is 



P d = l-exp(-N d ). 



(7) 



Qualitatively, Nd < 1 is required for a cluster to avoid a 
single, destructive encounter with a black hole. 

In the present work, our goal is to determine Pd using 
realistic calculations of cluster evolution: namely the Fokker- 
Planck solutions discussed above. As mentioned above, the 
method is linear and valid only for / < 0.15. Therefore, if 
we take fd — 0.15, then it is most appropriate to call Pd 
the probability for a strong collision. However, we also refer 
to Pd as the disruption probability according to its original 
definition. 

Since we are concerned with the probability that a clus- 
ter does or does not suffer one such collision, the range of 
possible final cluster states is broad: in other words, the evo- 
lution is stochastic. We therefore calculate collision proba- 
bilities using Monte Carlo simulations. 



Each simulation has 60 realizations. The collision his- 
tory for each realization is obtained by direct sampling of 
the relative velocity distribution and space density of black 
holes of mass Mbh within 10 initial tidal radii of the cluster. 
The disruption probability is calculated directly from the 
fraction of the 60 runs in which an encounter with / > 0.15 
occurs. Initially, each cluster has mass M c , King parameter 
Wo, galactocentric radius R and additional parameters (in- 
ternal mass spectrum: Murali & Weinberg 1997b,c) given in 
Table 0. 



3 CONSTRAINTS FROM GLOBULAR 
CLUSTERS 

3.1 Cluster evolution in a smooth halo:/bh = 

For fbh — 0, our calculations include only relaxation and 
post-collapse heating of the core. Table m gives the evapo- 
ration times tevap for clusters on circular orbits at R gc — 
16 kpc. The evaporation time is roughly 10 Gyr = tu for 
M c ~ 10 4 Mq, independent of concentration and roughly 
proportional to the initial mass of the cluster. The evap- 
oration times scale as t eva p(R) = t evap (R/lQ'kpc) for an 
isothermal halo. The lack of dependence on concentration is 
also seen in calculations presented by Lee et al (1991) and 
Gnedin & Ostriker (1997). 

^From these results, we conclude that clusters with ini- 
tial masses M c = 2.5 x 1O 4 M and 7.5 x 10 4 M Q provide ap- 
propriate specimens to study under bombardment by black 
holes: at this radius, they have t evap >> ts for fbh — 0; the 
low-mass cluster lets us probe the lowest values of Mhigh 
while the high-mass cluster is farther from the evaporation 
boundary, thus giving greater confidence in the conclusions. 



3.2 Collision probabilities in lumpy halos: fbh > 

Our main calculations depict the evolution of King model 
clusters on circular orbits at R gc = 16 kpc with a range of 
concentration and mass. We take halos with a range of fbh 
and Mbh = Mhigh, where Mhigh depends on the initial con- 
centration and mass of the cluster. Figures 1-3 compare the 
results of the Monte Carlo calculations with the analytic pre- 
dictions for the fixed cluster potential (equation ph. There 
appear to be statistical fluctuations in the Monte Carlo re- 
sults as well as a small systematic trend for 15% collision 
probabilities to lie below the fixed cluster predictions. Error 
bars are 68% confidence regions found using Bayesian meth- 
ods (as in Arras & Wasserman 1998) and are only statistical. 
Overall, however, the agreement is surprisingly good given 
that the analytic prediction neglects changes in the cluster 
potential due to internal evolution. 

Lower Pd might be expected in the simulations because 
two-body relaxation hardens the potential while mass loss 
reduces the cross-section for collisions. Although evolution 
should help clusters avoid strong collisions, Pd is only re- 
duced by roughly 10-20% for fbh ;$ 0.1. Agreement is best 
at the smallest and largest fbh- Agreement at fbh — > is 
trivial, as Pd — > in that limit. The agreement at larger 
fbh, where Pd — > 1, arises because the expected time to the 
first destructive encounter is small, so evolution (and the re- 
sulting increase in concentration) is relatively unimportant. 
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Table 1. Cluster Initial Conditions 



Structure 



Adopted values 



M c total mass 

Wq King concentration parameter 

R c cluster limiting radius 

t r h initial half-mass relaxation time 



2.5 x 10 4 ,7.5 x 1O 4 M 

Wo = 3,5,7 

set to tidal limit rt 

2 - 6 X 10 9 yr 



Internal mass spectrum 



(i mass spectral index: N(m) oc m^l 3 (i = 2.35 (Salpeter) 

mi lower mass limit mj = O.IMq 

m u upper mass limit m u = 2.0Mq 



Orbit: 



circular at 16 kpc 



Table 2. Evaporation times at 16 kpc (in Gyr) 



M C (M Q )/W 
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Figure 1. Probabilities for 15% encounters for Wq = 3 clus- 



ters with indicated masses at 16 kpc where M b 



M h 



igh 



Solid 



squares with associated error bars show the results of Monte Carlo 
simulations; dashed line shows the analytic prediction determined 
from equations (6) and (7). 



Figure 2. Probabilities for 15% encounters for Wo = 5 clus- 
ters with indicated masses at 16 kpc where M bh = M high . Solid 
squares with associated error bars show the results of Monte Carlo 
simulations; dashed line shows the analytic prediction determined 
from equations (6) and (7). 



Nevertheless, the agreement between the analytic formulae 
and simlulations is surprisingly good even for intermediate 
values, fbh ~ 0.05, where Pd ~ 0.3 — 0.5. 
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Figure 3. Probabilities for 15% encounters for Wo = 7 clus- 



ters with indicated masses at 16 kpc where Mf, 



M. 



high 



Solid 



squares with associated error bars show the results of Monte Carlo 
simulations; dashed line shows the analytic prediction determined 
from equations (6) and (7). 



Figure 4. Probabilities for 15% encounters for Wo = 5 clusters 
with indicated masses at 16 kpc where Mbh = 2M/ lig i l for Mhigh 
used in figures 1-3. Solid squares with assocated error bars show 
the results of Monte Carlo simulations; dashed line shows the 
analytic prediction determined from equations (6) and (7). 



3.2.1 Independence of Mbh > Mhigh 

To see how cluster evolution affects Pd at higher Mbh, we re- 
peat the above calculations for Wo = 5 and M^ = 2Mhi g h- 
Figure 4 shows the results. They agree well with the results 
for Wo — 5 presented in the previous section. We conclude 
that evolutionary effects do not destroy the high-mass scal- 
ing derived in the fixed cluster approximation. This, in fact, 
is not surprising since we have restricted our attention to 
clusters with t ev > in- 



3.2.2 Approximate behavior for Mbh < Mhigh 

In the above analysis, we have conservatively adopted a very 
large choice for Mhigh to ensure that our calculations obey 
the high-mass formalism in the strictest sense. However, Ar- 
ras and Wasserman (1998) have shown that the tidal limit 
formula for mass loss provides a good approximation even 
when bd < rt since clusters have fairly extended, loosely 
bound halos. Typically we expect the tidal formula to re- 
main fairly accurate even for impacts just outside the core 
radius (bd/ft ~ 0.1 for Wo = 5). This approximation is worst 
at low concentration because the mass distribution is more 
extended. The approximation improves as the cluster evolves 
and becomes more concentrated. For example, bd = 0.25r t 
encloses roughly 50% of the mass in the Wo — 3 cluster and 
roughly 80% of the mass in the Wo — 7 cluster. 

For fixed N d and V re i >max - 2.5{V re i), M bh oc b d . 
Thus the approximation significantly improves the limit- 
ing Mbh which our calculations explore. In particular, if 
we adopt bd — 0.25rt, then our calculations are valid for 
Mbh = 0.0625Mhi g h- In other words the high-mass scaling 
obtains for Mbh « Mhigh which we have defined above. 



Figure 5 shows the collision probabilities calculated using 
this approximation. The results agree well with the analytic 
predictions and indicate that the high-mass or destructive 
regime obtains down to Mbh ~ 10 6 Mp 



i @ . 



3.2.3 Radial scaling of collision probabilities 

The collision probabilities enumerated for 7? = 16 kpc can be 
scaled approximately to larger radius by keeping Nd fixed. 
Since Phaio oc R~ 2 and K oc R 2 from equation (El), then Nd is 
constant with radius for fbh oc R. The scaling is approximate 
because the intrinsic evolutionary rate of a cluster varies 
with radius due to the variation in dynamical time: tdyn °c R 
for a tidally limited cluster in an isothermal halo. However, 
the above calculations fall into the regime t evap > tu > T^ 1 , 
so the scaling is strong. Additional calculations verify the 
scaling. 



3.2.4 Eccentricity dependence 

Clusters on eccentric orbits in the Galaxy experience time 
variation in Td, the rate of destructive encounters, since the 
number density of black holes varies along their orbits. The 
galactic tidal force also varies along a cluster orbit; this 
can be accounted for roughly by assuming that the clus- 
ter is tidally limited at the pericenter of its orbit. The inte- 
grated number of destructive encounters can then be written 
Nd — Nd(e = 0)C(e), where e — 1 — J/J max (E) is a measure 
of the eccentricity of an orbit with energy E and angular 
momentum J, and Nd(e — 0) is the number of destructive 
encounters for a circular orbit with energy E (and angular 
momentum J max (£))- 
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Figure 5. Probabilities for 15% encounters for Wo = 5 clusters 
with indicated masses at 16kpc where M^h = 0.0625M/ t j 9 / t for 
M high used in figures 1-3. Solid squares with assocated error bars 
show the results of Monte Carlo simulations; dashed line shows 
the analytic prediction determined from equations (6) and (7). 



Figure 6. Mean final mass and concentrations of Wo = 3 clusters 
which do not suffer 15% collisions. Top row shows results for 
clusters with Mi„n = 2.5 X 10 4 Mq; bottom row shows results 
for clusters with M init = 7.5 X 1O 4 M . 



Whether or not Nd increases or decreases with increas- 
ing eccentricity depends on a competition between a higher 
number of encounters at smaller R, and more time spent 
at large R near apocenter. By numerical evaluation of C, 
we have found that C < 1, so that the penalty of spend- 
ing more time at apocenter with the accompanied small en- 
counter rate is the dominant factor. A convenient analytical 
formula for C(e) can by found by expanding the integral 
Nd{t) — J dt'Td(t') for small e and integrating t' over an 
integral number of orbital periods - a good approximation 
if many orbits have been traversed. Under these conditions 
we find 



C(e) 



1 



? e V2 + 33 ( 3 /a 

2 24 v ; 



(8) 



The probability of survival is larger for clusters on ec- 
centric orbits than for clusters on circular orbits, given E. 
For example, consider two clusters with the same orbital en- 
ergy, one on a circular orbit at R = 16 kpc and a second on 
an eccentric orbit with pericenter at R p — 8 kpc; in this case 
Nd is smaller for the eccentric orbit by about a factor of two. 
In this case, the cluster on the eccentric orbit has a higher 
rate of internal evolution because of its smaller pericenter, 
an effect which must also be taken into account; we discuss 
qualitatively in §4. Finally, we note that the average value of 
the correction factor for an isotropic distribution of angular 
momenta is (C(e)) ~ 0.5 for a fairly large range of orbital 
energies. 



3.3 Properties of evolved clusters 

We examine the basic properties of the clusters which do not 
undergo 15% collisions. Thus we describe the effect of the 



weaker encounters on cluster evolution. Ideally we would like 
to determine the Green's function or probability amplitude 
for evolution from a given initial state to a given final state 
(e.g. Arras & Wasserman 1998). Of course, with 60 realiza- 
tions per run and only the fraction 1 — Pd(fbh) not suffering 
strong collisions, we can only understand the distribution of 
final states very approximately. To do so, we simply examine 
the mean mass and concentration of the 'surviving' clusters. 

Figures 6-8 show the mean mass and concentration for 
clusters in the runs discussed above. In each case, remaining 
mass decreases montonically with fbh- As mentioned above, 
the black hole flux leads to both weak and strong encounters: 
clusters which do not suffer strong collisions still lose mass 
through weak encounters, so the final mass will be less than 
that of an isolated cluster. 

The evolution of concentration c = log rt/r c behaves 
somewhat differently in each case. At low mass, the relax- 
ation rate is enchanced by weak encounters for all initial 
concentrations. However, the Wo = 3 clusters have not yet 
entered core collapse so c increases montonically with fbh- 
For Wo — 5 and Wo — 7, c decreases with fbh because 
the core reaches the post-collapse stage of expansion more 
rapidly due to the external heating. 

At high mass, the effect differs because of the longer 
intrinsic relaxation time. For initial Wo — 3, c decreases 
with fbh because heating has little effect on core evolution, 
tending only to decrease rt. For Wo — 5, c increases with fbh 
due to the acceleration of core evolution by external heating. 
For Wo = 7, external heating accelerates core evolution past 
collapse and into expansion, so that c decreases with fbh- 
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Figure 7. Mean final mass and concentrations of Wo = 5 clusters 
which do not suffer 15% collisions. Top row shows results for 
clusters with Mi n u = 2.5 X 10 4 Mq; bottom row shows results 
for clusters with M init = 7.5 X 10 4 M Q . 
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Figure 8. Mean final mass and concentrations of Wo = 7 clusters 
which do not suffer 15% collisions. Top row shows results for 
clusters with Mi„u = 2.5 X 1O 4 M0; bottom row shows results 
for clusters with Mi n u = 7.5 x 1O 4 M0. 



4 DISCUSSION 

We have examined in detail how typical globular clusters 
evolve in a halo which contains a population of massive 
black holes with Mbh ~ 1O 6 M0. Our main goal was to estab- 
lish probabilities for strong collisions between clusters and 
black holes in the high-mass limit using Fokker-Planck cal- 



culations in order to combine effects of internal relaxation, 
binary heating and black hole shocking. 

Our results show that evolution does not radically alter 
collision probabilities determined in the fixed cluster approx- 
imation (equation 0). This result is not surprising given our 
approach: we first determine initial conditions which do not 
lead to significant evolution in the absence of black holes; we 
then include black holes of some mass Mbh and abundance 
fbh in calculations with these initial conditions. Evolution 
can help reduce Pd, but weaker encounters tend to acceler- 
ate core collapse and evaporation by removing mass from 
the halo. The differences that do arise are relatively small 
and therefore will not affect any conclusions we draw below. 

In calculating Pd, we have only considered clusters on 
circular orbits. As discussed in §3.2.4, clusters on eccentric 
orbits have lower collision probabilities. However, decreasing 
the pericenter also shortens the evaporation time for a clus- 
ter in isolation. For example, Table 2 shows that a cluster 
with mass M c = 2.5x 10 4 Mq has an evaporation timescale of 
about 20 Gyr assuming a circular orbit of radius R — 16 kpc; 
if its pericenter were R p = 8 kpc, its evaporation time would 
be about 10 Gyr. Moreover, disk shocking becomes effective 
for pericenter radii R p < 8 kpc, further reducing the odds of 
survival for M c = 2.5 x 10 4 M Q even if fbh — 0. To be con- 
servative, we could avoid these complications by focusing on 
clusters of larger mass, resulting in correspondingly larger 
values of the minimum black hole mass on which we can 
place constraints in the high-mass limit: for a fixed evapora- 
tion time, and assuming tidally limited clusters, we should 
scale M c oc Rp 1 and Mhi g h oc M c R p oc Rp , approxi- 
mately. If we choose to be bolder, we could keep M c fixed, 
restrict ourselves to acceptable values of R p (e.g. R p > 8 kpc 
for M c = 2.5 x 1O 4 M ) and rescale M high oc R^ /3 . Either 
way, the limits on fbh are unlikely to change by more than 
a factor of two (see §3.2.4). Because we have only rigorously 
considered clusters on circular orbits, we will address any 
remaining ambiguities in defining Mhigh in our next paper 
(Arras et al 1999), where we shall consider the evolution of a 
realistic population of globular clusters and work our way up 
from the regime of low black hole masses, thereby obtaining 
limiting values of fbh as a function of Mbh', where that curve 
asymptotes to the Mth-independent bound found here will 
determine Mhigh- 

We have calculated Pd assuming fd = 0.15; this restric- 
tion was imposed by the limitations of our linear pertur- 
bation procedure. We can scale the results to higher fd to 
consider collisions more likely to disrupt the cluster, i.e., 
fd ~ 0.5. Here nonlinear effects become important but the 
scaling approximately holds. In linear theory, N d oc fbh/^ffd 
from equation (6) while M high oc y/Ja- Therefore, for fixed 
fbh, Nd oc \j\fjd and the values of Pd calculated above de- 
crease correspondingly for fd > 0.15. For fixed Nd and Pd, 
fbh oc ■s/Jbh - (Note that if N(f' d ) ~ 1, then for f d < f' d , 
Nd(fd) ~ \/f'd/f d )> w hich would only result in a fractional 
mass loss ~ yfdfd < f'd': large fractional mass loss in a 
single encounter is always more likely than in numerous en- 
counters each with smaller fractional mass loss.) 

We can also scale to different galactocentric radii us- 
ing the approximate relationships fbh oc R (see §3.2.3) and 
Afhigh oc fl /2 Mc /3 R 1/3 , which becomes M high oc f^ 2 R~ iri 
for fixed evaporation time (see §3.2.4); including disk shock- 
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ing - which we have neglected here - ought to raise the 
value of Mhigh somewhat because only higher mass clusters 
survive. The constraints on fbh become stronger for smaller 
galactocentric radii, but we expect f b h to be nonuniform as 
a consequence of dynamical friction, so limits at relatively 
large R are easiest to interpret. 



4.1 Interpreting the collision probabilities 

Our calculations in smooth halos, fbh = 0, combined with 
the observational picture serve as a guide for interpreting 
the collision probabilities. Observationally, the similarity of 
globular cluster luminosity functions over a range of galaxy 
environments (e.g. Harris 1991) may reflect the formation 
process and suggests that cluster populations are relatively 
unevolved for M c ~ 1O 5 M0. Moreover, the flatness of the 
distribution of globular cluster luminosities, dN/dL c - and 
hence dN/dM c - at low L c in the Milky Way (e.g. Ashman & 
Zepf 1998) suggests that globular clusters with M < 10 5 M Q 
are not whittled away rapidly, but are relatively long-lived 
since, otherwise, the rate of destruction would far exceed the 
rate of production and the distribution would fall off drasti- 
cally. Our calculations with fbh = appear to be reasonably 
consistent with this expectation. In particular, these calcu- 
lations show that clusters on circular orbits at 16 kpc survive 
beyond a Hubble time for M > 1O 4 M0, incurring roughly 
25% mass loss for M c = 2.5 x 1O 4 M0 and roughly 5% mass 
loss for M c = 7.5 x 1O 4 M0. A stringent view, therefore, re- 
quires that halo black holes leave these clusters relatively 
unscathed. 

Our numerical simulations show that the probability 
of a collision with fd — 0.15 is about 30-40% for fbh ~ 
0.025, about 50% for f bh ~ 0.05, and about 80% for f bh » 
0.1; assuming that the collision probability in this range of 
fbh is excessive, we can rule out Mbh ^ 1-3 x 1O 6 M0 for 
M c = 2.5 x 10 4 A/ Q and M bh > 2.7 x 10 6 Af Q for M c = 
7.5 x 10 4 M Q for f bh > 0.1. From the scaling N d oc fbh/VTd, 
equal collision probabilities for fd = 0.5 imply the values 
fbh = 0.05, fbh = 0.09 and fbh = 0.18, respectively; since 
Mhigh oc t/J2 (see Section 2.2, especially eq. [3]), these limits 
apply to Mbh > 2.4 x 10 6 M for M c = 2.5 x 10 4 M and 
M bh > 4.9 x 10 6 M Q for M c = 7.5 x 1O 4 M . These are 
somewhat higher but still very restrictive. 

Black hole collisions do not only destroy clusters, but 
also whittle away their masses. For low fdNd, the fractional 
mass loss endured by surviving clusters is approximately 
fdNd in the tidal limit (Figs. 6-8 and Arras & Wasserman 
1999). Thus, we expect that for large fdNd, the mean mass 
per cluster declines like exp(— fdNd), and the total mass 
in clusters declines like exp[— (1 + fd)Nd] when cluster de- 
struction is taken into account, in the tidal limit. To a first 
approximation, the evolution of the distribution of clusters 
must account for a steady advection downward in mass as 
well as destruction. If N(M, t)dM is the number of clusters 
with masses between M and M + dM at time t, then 

dN ^^-YdN { M,t )+ Ydfd d[N ^! )M] (9) 



i)t 



dM 



represents a simple advection-destruction model appropriate 
for the tidal limit. The solution to this equation is 

N{M,t) = exp[-r d (l - f d )t]N(Mexp(f d T d t),0). (10) 



According to this solution, the total number of clusters de- 
clines oc exp(— Tdt) — exp(-Nd) and the total mass re- 
maining in clusters after time t is oc exp[— F d (l + fd)t] = 
exp[-(l + f d )N d }. 

For a given black hole mass Mbh, the tidal approxima- 
tion holds up to some maximum cluster mass, Mt ; if atten- 
tion is restricted to the tidal approximation, which considers 
only collisions with impact parameters outside n, this max- 
imum mass is M t ~ (a few) x 10 4 M Q for M bh = 10 6 M Q . 

Black hole collisions still destroy and whittle away clus- 
ters with masses above M t . Consequently, clusters destroyed 
and chiselled away at masses below Mt are replaced, to a 
degree, by clusters originally at masses above Mt- Our cal- 
culations in the tidal limit cannot describe this evolutionary 
process; to do so requires calculations of what happens as a 
consequence of black hole collisions at b < r t , which we have 
excluded in this paper (but are in the midst of calculating, 
and will publish separately). However, from earlier work (e.g. 
Arras & Wasserman 1999), we already know that the tidal 
approximation continues to describe the mean mass loss by 
a cluster to within 20-30% for somewhat smaller impact pa- 
rameters, b > (0.1 — 0.2)7-4. A quantitatively reasonable in- 
terpolation formula for the fractional mass loss of a cluster 
due to a collision with a black hole passing within impact 
parameter 6 of the center of a cluster at relative speed v re i 
is 



f(b,v rel ) = 



f(0,Vc)(Vc/Vrel) 2 

[1 + {b/boYf 



(11) 



where /(0, V c ) is the fractional mass loss at zero impact pa- 
rameter for relative velocity v re i = V c ; for a tidally limited 
cluster at galactocentric radius r, ,f(0, V c ) oc M b 2 h /M 4/3 r 2/3 . 
The numerical value of /(0, V c ) can be found using the re- 
sults of Arras & Wasserman 1999 (see fig. 3). The value of 60 
is fixed by requiring the rate of destructive encounters found 
using eq.hll give the correct tidal limit. 

For a given black hole mass, there is a new critical 
mass Md above which clusters become considerably more 
immune to destruction. Typical penetrating encounters are 
nondestructive when the cluster mass is large enough that 
/(0, V c ) < fd- Using eq. (Ill]), we estimate 



M d ^2x 10 5 M Q ( Mbh 



V 10 6 M 



\',/2 



8 kpc 
R 



1/2 



r*- 



The evolution of clusters with M > Md in the face of black 
hole collisions is described poorly by the tidal approximation 
for two reasons. First, the rate of destructive encounters is 
low, because /(0, V c ) drops below fd for M > Md- Second, 
the rate at which cluster advect downward in mass is slowed 
because of the same cutoff in fractional mass loss. We can 
generalize eq. fl9J) to account for the different behavior at 
masses above and below Md', the resulting equation is 



dN{M,t) 
di 



+ IV, 



T d I(M/M d )N(M,t) 

d[N{M, t)MH(M/M d )] 



dM 



(13) 



where, from the scalings implied by eq. (IllT) , we infer that 
I(z) -* 1 for 2 < 1 and I{z) -» for z > 1, and H{z) -> 1 
for 2 < 1 and H(z) ~ z~ 2/3 for 2 > 1. (The exact M/M d 
dependences can be computed from eq. [[11| given the distri- 
bution of V .) 
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The tidal limit is recovered for essentially the entire 
range of globular cluster masses for sufficiently large M b h- 
The most massive Galactic globular cluster, lo Cen, has 
M c w 2.4 x 10 6 M (e.g. Ashman & Zepf 1998, Mandu- 
shev et al. 1991), and M d > 2.4 x 10 6 M for M bh > 
5.2 x 10 6 M Q (7?/8kpc) 1/3 (0.5// d ) 1/2 from eq. @. Since 
cluster masses decrease with time as a result of black hole 
perturbations, Mbh must be somewhat larger still for the 
tidal approximation to hold for all time. 

According to the advection-destruction solution, eq. 
(P_0|), the cluster mass distribution evolves self-similarly with 
time in the tidal limit, which means that the initial condi- 
tions must have resembled the distribution seen today, ex- 
cept shifted to larger mass. Thus, we cannot use the shape 
of the distribution to derive a bound on fbh, without addi- 
tional information on the initial mass function for clusters, 
and without including other processes, such as evaporation 
and (at R < 8 kpc) disk shocking, rl Requiring that the total 
mass lost by clusters not exceed about 100 times the total 
mass presently contained in them, so the halo star popula- 
tion is not primarily due to ejecta from clusters (e.g. Klessen 
& Burkert 1996 and references therein), implies Nd J$ 3.1, 
or fbh % 0.36(-R/16 kpc), for f d = 0.5. This limit is fairly 
conservative, since metallicity differences suggest that the 
bulk of halo field stars did not originate from globular clus- 
ters (Harris 1991). Moreover, clusters on nearly circular or- 
bits are somewhat likelier to be disrupted than clusters on 
elongated orbits, so we should expect ejected stars to have 
a tangentially-biased velocity ellipsoid in the inner Galaxy, 
which is not observed for the halo stars (Beers & Sommer- 
Larson 1995). 

For smaller M b h, for which Md falls in the range of 
present (and past) cluster masses, we can use eq. (|13j) to get 
a rough idea of how the cluster distribution function evolves 
as a consequence of bombardment by massive black holes. 
(We shall present more realistic and accurate models, as well 
as statistical analyses, in a subsequent paper.) Qualitatively, 
clusters with masses below Md ought to lose mass and be 
destroyed rapidly. Clusters with masses above Md are less 
prone to destruction, and lose mass more slowly. Although 
the advection of clusters from masses above M d to below 
Md ought to replenish the supply of low mass clusters lost 
to destructive collisions, the slowness of the process results 
in an overall truncation of N(M, t) at low values of M. 

To explore the stability of the presently-observed mass 
distribution of clusters, we used eq. ( |l3[ ) to compute the 
evolution of a distribution that is originally 



N(M,0) = 



Mo(l + M/M ) 



(14) 



with Mo — 3 x 10° Mo, up to a maximum mass M max = 
1O 7 M0; the calculations assumed fd — 0.5 and Md = 
6 x 10 5 M Q (corresponding to M bh = 2 x 10 6 M ). The orig- 
inal cluster mass distribution given by eq. (Il4j) is flat below 



Since the evaporation time scales oc M, the rate of evapo- 
rative mass loss, M ev , is approximately independent of mass, 
and the tidal limit solution with evaporation is N(M, t) = 
cxp[-(l - f d )t]N(Mi(t),0), with Mi(t) = Mexp(f d r d t) + 
(M.„ IT j f j)\exp(f jT j£\ — 1]. Evaporation can be included sim- 
ilarly in solving eq. (J13|) . In deriving limits on fbh, we neglect 
evaporation, which would make our constraints slightly tighter. 



Mo and oc M~ above, in reasonable agreement with the ob- 
served luminosity function for globular clusters in the Milky 
Way (and other galaxies), assuming a constant M/L (a 3 
(e.g. Ashman & Zepf 1998). The results for MN(M,t), the 
distribution of clusters in In M, is shown in Fig. 9a; the loga- 
rithmic slope d In N(M, t)/d In M is shown in Fig. 9b. These 
figures illustrate the truncation of the distribution at low 
masses, and show how the shape of the distribution tends 
to evolve with time. From Fig. 9,b, we see that the shape 
begins to deviate significantly from the original one after 
Tdt — Nd ~ 2 or 3. Although this result is preliminary, we 
expect that more sophisticated analysis will lead to a similar 
bound. 

We can also examine the evolution of the globular 
cluster system from presumed initial conditions using eq. 
(I13j). As in the tidal case, we remark that the evolution- 
ary sequences computed here only include the effects of 
black hole perturbations, not the better established effects 
known to operate, such as evaporation and disk shocking. 
For illustrative purposes, we adopted N(M, 0) oc M~ 2 for 
10 3 M Q < M < 10 7 M Q (e.g. McLaughlin 1999). The results 
are shown in Fig. 9c, d. As can be seen, N(M, t) evolves to 
a form reminiscent of today's cluster mass distribution over 
a time span Tdt ~ 10, which corresponds to w 10 Gyr for 
fbh ~ 1.1(-R/16 kpc). However, it is also clear from this fig- 
ure, and from Fig. 9a, b, that the shape of N(M, t) evolves 
somewhat more rapidly, on a timespan TdAt ~ 2 or 3, so 
the observed distribution is not long-lived. We note that the 
total mass in clusters drops by about a factor of ten between 
Tdt — and Tdt — 10 for this model, which, although sub- 
stantial, does not violate any observational constraints (cf. 
Klessen & Burkert 1996). 

From these preliminary investigations, we conclude that 
while it is possible to find initial conditions N(M, 0) that 
evolve to the observed distribution of cluster masses in ~ 10 
Gyr even for fbh ~ 1, the evolutionary timescale is rather 
short, so the distribution is rather unstable. Although a 
more precise treatment of the bounds on fbh that can be 
derived from requiring stability of the cluster mass distri- 
bution needs to be done (and will be presented by us else- 
where), we are confident that the observed distribution is 
unstable unless Tdt — Nd J$ 3 everywhere in the halo. 
Adopting this rather conservative limit, we conclude provi- 
sionally that stability of the cluster mass function requires 
fbh <0.34(i?/16kpc). 

We therefore propose two different bounds on fbh- The 
less restrictive of the two is the stability bound derived 
above. For direct comparison with the disk heating bound, 
we shall rescale all constraints to R — 8 kpc; then stability 
implies fbh ^ 0.17(i?/8kpc). This bound applies, strictly 
speaking, to Mbh = 2 x 10 6 M© . At much larger black hole 
masses, where the tidal approximation is valid, we found a 
similar bound based on mass loss from clusters to the halo, 
and we adopt f bh < 0.17(R/8 kpc) for all M bh > 2 x 10 6 M Q . 
The more restrictive bound is obtained by assuming that a 
hypothetical cluster with M c — 2.5 x 10 4 Mq on a circular or- 
bit at 16 kpc should have no worse than a 50% chance of sur- 
vival. Using the scalings given above to compute bounds at 
R = 8 kpc, we find f bh < 0.025 for f d = 0.15 and f bh < 0.05 
for f d = 0.5. If we keep M c = 2.5 x 10 4 M Q fixed, then these 
limits apply to M bh > 1 x 10 6 M© and M bh > 2 x 10 6 M Q , re- 
spectively. Since the evaporation time falls to about 10 Gyr 
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Figure 9. Panels (a-d) show the cluster population mass function 
and its power law index for a range of times. The initial condition 
N(m,0) = (1 + M/M )- 2 /M was used in figures (a) and (b), 
and the initial condition N(m, 0) = 1/M 2 was used in figures (c) 
and (d). The eleven curves in each figure correspond to the times 
T^t = Nd = 0, 1, ..., 10. In figures (a) and (c), time increases as 
the curves go from top to bottom, while in figures (b) and (d) the 
opposite is true. The sharp edges seen at large masses occur where 
the mass function is zero in the advection/destruction model. 



at R = 8 kpc for M c = 2.5 x 1O 4 M and f bh = 0, we might 
prefer to keep the evaporation timescale fixed at 20 Gyr, in 
which case our derived limits on f b h apply to slightly larger 
black hole masses, M bh > 2 x 1O 6 M and M bh > 3 x 1O 6 M 
for fd = 0.15 and fa — 0.5, respectively. For our conserva- 
tive bounds, we adopt f bh < 0.05 for M bh > 3 x 1O 6 M ; 
for our most stringent bounds, we adopt f b h < 0.025 for 
M b h ^lx 1O 6 M . Figure 10 shows these limits along with 
the upper limit 



M bh < 4.4 x IO 7 M 



8 



R 



In A / V 8 kpc 



(15) 



for In A = 8 imposed by requiring that dynamical friction 
be incapable of dragging black holes inward in 10 Gyr (see 
equation [7-27] in Binney & Tremaine 1987), and the disk 
heating constraint f b uM b h <2x IO 6 M (Lacey & Ostriker 
1985). The figure shows that the globular cluster constraint 
forbids considerable portions of M b h — ft>h space allowed by 
disk heating. 
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Figure 10. The constraints on the mass and fraction of black 
holes at 8 kpc in the halo. The solid horizontal line shows the up- 
per limit on logMj,^ from dynamical friction. For logM bh > 7.7, 
black holes spiral into the Galactic center. The solid diagonal 
line shows the disk heating constraint: log /j^ Mj,^ > 6.3 over- 
heats the disk. The solid box labeled stability shows the con- 
straint obtained from the stability requirement for the globular 
cluster population, log f b h < —0.77 for log M b ^ > 6.3. The dashed 
box shows the more conservative constraint for the weak cluster 
evolution hypothesis log f b h < —1.3 for logM^ > 6.3 and the 
solid box labeled weak evolution shows the most stringent bound, 
log/fch < —1-6 for logMj,/j > 6. The shading indicates regions 
which arc disfavored by the combination of constraints. 



APPENDIX A: SECOND-ORDER CHANGE IN 
DF IN IMPULSIVE ENCOUNTER 

We derive the second-order change in the cluster DF due 
to an impulsive, tidal encounter with a perturber. As men- 
tioned above, the resulting expression is equivalent to a 
Fokker-Planck equation for the change in the DF and there- 
fore consists of an advection term and a diffusion term. Our 
derivation reveals an error in previous treatments (Kundic 
& Ostriker 1995; Gnedin & Ostriker 1997; Murali & Wein- 
berg 1997a) which results from improper integration over 
velocity coordinates in projecting the advection and diffu- 
sion coefficients to energy space. As shown in the derivation 
and discussion below, the problem arises because one can- 
not integrate over the entire range of angular coordinates in 
velocity space: failure to restrict the integration domain cor- 
responds to including fictitious transitions to bound states 
from unbound states that are, in reality, unoccupied initially. 
Mathematically, there is a ^-function implicit in the kinetic 
equation which has previously been ignored. 

For any impulsive enounter with a perturber, the equa- 
tion, 

/ne»(r',v') = /drdv5(r-r')<5(v-v' + Av(r))/(r,v) 
= /(r',v' -Av(r')); (Al) 
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gives an exact expression for the new DF, / neTO , in terms of 
the old DF /. The 5-function in time indicates that the per- 
turbation is impulsive; the (^-function in position indicates 
that particles do not move during the perturbation; and the 
5-function in velocity defines the position-dependent velocity 
impulse (with respect to center of mass), Av(r), imparted 
to a particle by the perturbation. Consequently, the new 
DF is the old DF with velocity bins shifted according to the 
position-dependent velocity impulse. 
The resulting mass loss 



S.\J = I dr'dv' f new (r',v'), 

'|v'|>l! c (l-') 



(A2) 



is the integral over all particles whose new velocity is gre ater 
than the escape velocity v e (r). Substituting equation ( |Al[ ) 
for /„ 



5M ■ 



we can show that 



dr / dv9(|v + Av(r)| - v e (r))f(r,v), 



(A3) 



where we have dropped the primes on the spatial coordi- 
nates. This is precisely equation (39) of Chernoff, Kochanek 
& Shapiro (1987) (hereafter CKS); therefore the present 
treatment of individual collisions is equivalent to that used 
by Arras & Wasserman (1998). 

To use a one-dimensional, phase-space method like the 
standard Fokker-Planck calculation, we must project the 
new DF into energy space by integrating over all other co- 
ordinates. The projection 

lQiv 2 P(E)f new (E) = / drdSl-u V / 2( J B - $)/ new (r, v) = 

/ drdn v y / 2(E - $)/(r, v - Av(r)) (A4) 

J v — Av|<u e (r) 

where the phase-space volume 



P{E) = / drr 2 v / 2[E-${r)]. 



(A5) 



The second integral in equation (A.4) has a limited range 
of integration because the perturbation can transport un- 
occupied, unbound states to bound states: areas where the 
unperturbed DF vanishes, i.e. E > Et, define the excluded 
regions. To reiterate, previous treatments have overlooked 
this subtlety in deriving the kinetic equation, thereby ob- 
taining a factor of 2 overestimate in the mass loss. 

To derive the second-order change, we expand f n ew in 
a Taylor series about the unperturbed DF: 

f new (r, v) = /(r, v - Av(r)) a /(r, v) - |£ ■ Av(r) 

4 Av(r) • £L ■ Av(r) - (A6) 

This has the standard form of the velocity-space Fokker- 
Planck equation. Substitution into equation (A.4) yields an 
equation for the change in the DF, Sf. 

We will now briefly outline the calculation of Sf. Let 
the direction of the relative velocity be along the z— axis 
and the position of a star in the cluster given by spherical 
radius r and the cosine of the angle with respect to the z 
axis, /i = cos 6. The magnitude of the tidal velocity kick 
with respect to the center of mass of the cluster is then [see, 



e.g., Arras & Wasserman (1998)] 

. . . 2GM bh ry/l-^ 
Avr) = rrrr • A7 

A V re l 

The region of phase space for which |v — Av| < v e (r) has 
been discussed in both CKS and Arras & Wasserman (1998). 
They find that, depending on the energy E and the position 
r and fi, the velocity angle cosine /j, v = v- Av/wAw may only 
extend over a restricted interval instead of (—1, 1). Hence, 
spherical symmetry is broken and the first order term in Av 
in equation A6 no longer integrates to zero. The end result 
is that there will be two different mathematical forms for Sf 
depending on the energy E. Define the critical energy 

2G MbhT peak 



E, 



— E n 



~Ve ^fpeak ) 



(A8) 



b 2 V rel 

where r pca k is the radius at which rv e (r) reaches a maxi- 
mum. For energies E < _E C rit, the velocity cosine fi v is in 
(—1,1) for all values of r and (J, and one would obtain the 
"standard" results for Sf. For E > E CI n, on the other hand, 
fl v runs over a restricted range and a different mathematical 
expression for Sf is obtained. 

It is convenient to write Sf in the quasilinear form 



J/(B)=[lftr>P(B)]- 1 ^ 

For energies E < £ C rit, we find 

64^ 2 df{E) ( GM bh \ 



(A9) 



9 dE \ b 2 V rA J 

<p- x (E) 

drr 4 (2[E ~ <j>{r)]f /2 , 



(A10) 



where <f>~ 1 (E) is the apocenter of a radial orbit with energy 
E. For E > -E^it, a more complicated expression results, 
given by 

64ir 2 df(E ma ^ '™-' 



F(E) = +- 



I) 



dE 



( GMbh Y 



drr v e (r) 



-8tt 



2 djyEmax 



dE 



AE) 



drr v e {r) 



,(s) 



"i \J-Jmax 

1 2GM bh r 



b 2 V rel 
1 1 b 2 V, 



E) 2 iio(r,E) 



-l V 



(E max - E) || - 0o (r, E) + i sin(2e„ (r, E) X 
«r : ur2GU b tr {Ema *- E)3 (h 6 ^v) 



+ 



1 a ,( 2GM bh r \ 2 f 1 3 \ 



(All) 



Here E — v 2 /2 + cj>(r) is the energy (per unit mass), E max — 
4>(rt) and E m in = 0(0), r± is the tidal radius of the cluster, 
/ is the pre-collision distibution function, r max and r m in are 
defined by the equation 



\1-Jmax -L'J — ^e\rmin,max jrmin.max : \ J 

cos(8o(r,E)) is defined by the 



2GM bh 

and the cosine fio(r,E) 
equation 



1/2 



" \-E-Jmax ^ } • 



(A13) 



2GM bh rv e (r) 
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For most of the range of E, the above equation gives the 
same result as if there were no restrictions on phase space. 
Only for energies within roughly Av/u of Emax do the re- 
strictions make a difference. However, the flux at the bound- 
ary, which is the change in the mass of the cluster is signif- 
icantly altered; we find 



F{E„ 



AM 



32tt 2 df(E„ 







[Emax) ( GM bh \ 

dE \bW re J 



drr v e (r) 



This is exactly a factor of two smaller than the answer ob- 
tained from ignoring the phase space restrictions. 

For reference, we give the value of K used in equation 
(2): 



K 



32*-" 2 1 
~ G M 



df(E„ 



dE 



drr v e (r). 



(A15) 



The above procedure is inherently linear and, therefore, 
has only a limited range of validity. The method fails because 
the new DF, /o + A/, becomes negative as we increase the 
mass loss. We find that negative excursions in the new DF 
become unacceptably large for dM/M > 0.15. 
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